
include("../src/qmcmary.jl")
using ..qmcmary

using LinearAlgebra

function Zonesite_tr(Nt)
    U = 1.5
    mu1 = 1.0
    mu2 = -1.0
    dt = 0.1*Nt
    Z = exp(-dt*U*0) + exp(-dt*U*0+0.5*dt*U-dt*mu1) +
    exp(-dt*U*0+0.5*dt*U-dt*mu2) + exp(-dt*U*1+0.5*dt*U*2-dt*(mu1+mu2))
    return Z
end


function Zonesite_hs(Nt)
    theta = 0.5pi
    phi = 0.5pi
    nx, ny, nz = cos(theta)*sin(phi), sin(theta)*sin(phi), cos(phi)
    U = 1.5
    mu1 = 1.0
    mu2 = -1.0
    #dt = 0.1
    #lam = acosh(exp(U*dt/2))
    #
    hslen = Nt
    hspool = ["-1", "+1"]
    for idx in Base.OneTo(hslen-1)
        hspool1 = hspool
        hspool2 = copy(hspool)
        for idx in Base.OneTo(length(hspool1))
            hspool1[idx] = hspool1[idx]*"-1"
        end
        for idx in Base.OneTo(length(hspool2))
            hspool2[idx] = hspool2[idx]*"+1"
        end
        hspool = vcat(hspool1, hspool2)
    end
    #
    Zp1 = 0.0
    println(hspool)
    for hscfg in hspool
        bseq = I(2)
        #println(hscfg)
        cfg = [parse(Float64, hscfg[(2*bi-1):2*bi]) for bi in Base.OneTo(Nt)]
        #println(cfg)
        for bi in Base.OneTo(Nt)
            c = cfg[bi]
            if (bi % 2) == 1
                #println(bi, " 1")
                dt = 0.1+0.0im
                lam = acosh(exp(U*dt/2))
            else
                #println(bi, " 2")
                dt = 0.1-0.0im
                lam = acosh(exp(U*dt/2))
            end
            clam = c*lam
            #M11 = e^{lam}(lam+lam*nz)/2/lam - e^{-lam}(lam*nz-lam)/2/lam
            M11 = exp(clam)*(1+nz)/2 - exp(-clam)*(nz-1)/2
            #M12 = (lam*nx-lam*ny*im)*e^{lam}/2/lam - (lam*nx-lam*ny*im)*e^{-lam}/2/lam
            M12 = (nx-ny*im)*exp(clam)/2 - (nx-ny*im)*exp(-clam)/2
            #M21 = (lam*nx+lam*ny*im)*e^{lam}/2/lam - (lam*nx+lam*ny*im)*e^{-lam}/2/lam
            M21 = adjoint(M12)
            #M22 = e^{lam}(lam-lam*nz)/2/lam - e^{-lam}(-lam-lam*nz)/2/lam
            M22 = exp(clam)*(1-nz)/2 - exp(-clam)*(-1-nz)/2
            cfg1 = [M11 M12; M21 M22]
            if (bi % 2) == 1
                dt = 0.1+0.1im
            else
                dt = 0.1-0.1im
            end
            tk = [exp(-dt*mu1) 0; 0 exp(-dt*mu2)]
            bseq = cfg1*tk*bseq
        #println(cfg1)
        #Zp1 += 0.5*det(I(2) + cfg1)
        end
        Zp1 += (0.5^Nt)*det(I(2) + bseq)
    end
    return Zp1
end


function Zonesite_hs_Gauge1(Nt)
    U = 1.5
    dt = 0.1
    mu1 = 1.0
    mu2 = -1.0
    ra = 2acosh(exp(U*dt/2))/U/dt - 3*rand()
    ia = 3*rand()
    @show ra ia
    #println(ra, " ", ia)
    rb = 0.75
    ib = 0.25
    #
    hslen = Nt
    hspool = ["-1", "+1"]
    for idx in Base.OneTo(hslen-1)
        hspool1 = hspool
        hspool2 = copy(hspool)
        for idx in Base.OneTo(length(hspool1))
            hspool1[idx] = hspool1[idx]*"-1"
        end
        for idx in Base.OneTo(length(hspool2))
            hspool2[idx] = hspool2[idx]*"+1"
        end
        hspool = vcat(hspool1, hspool2)
    end
    #
    Zp1 = 0.0
    println(hspool)
    for hscfg in hspool
        bseq = I(2)
        #println(hscfg)
        cfg = [parse(Float64, hscfg[(2*bi-1):2*bi]) for bi in Base.OneTo(Nt)]
        #println(cfg)
        wgt = 1.0
        for bi in Base.OneTo(Nt)
            c = cfg[bi]
            dt= 0.1
            y1 = exp(dt*U/2)
            y1a::ComplexF64 = y1^(ra + ia*im)
            pval = rb + ib*im
            M11 = 0.5*(c+1)*(y1 - pval*y1a) / (1-pval) - 0.5*(c-1)*y1a
            M12 = 0.0
            M21 = 0.0
            M22 = 0.5*(c+1)*(y1a*y1 - 1)/(y1a - y1) -
            0.5*(c-1)*(1 - pval - y1^2 + pval*y1a*y1)/(pval*y1a - pval*y1)
            #
            cfg1 = [M11 M12; M21 M22]
            tk = [exp(-dt*mu1) 0; 0 exp(-dt*mu2)]
            bseq = cfg1*tk*bseq
            wgt = wgt*(0.5*(c+1)*(1-pval) - 0.5*(c-1)*pval)
        #println(cfg1)
        #Zp1 += 0.5*det(I(2) + cfg1)
        end
        Zp1 += wgt*det(I(2) + bseq)
    end
    return Zp1
end


function Zonesite_hs_Quad2(Nt, phi)
    U = 1.5
    dt = 0.1
    mu1 = 1.0 - U*phi
    mu2 = -1.0 - U*phi
    #
    hslen = Nt
    hspool = ["-2", "-1", "+1", "+2"]
    for idx in Base.OneTo(hslen-1)
        hspool1 = hspool
        hspool2 = copy(hspool)
        hspool3 = copy(hspool)
        hspool4 = copy(hspool)
        for idx in Base.OneTo(length(hspool1))
            hspool1[idx] = hspool1[idx]*"-1"
        end
        for idx in Base.OneTo(length(hspool2))
            hspool2[idx] = hspool2[idx]*"+1"
        end
        for idx in Base.OneTo(length(hspool3))
            hspool3[idx] = hspool3[idx]*"-2"
        end
        for idx in Base.OneTo(length(hspool4))
            hspool4[idx] = hspool4[idx]*"+2"
        end
        hspool = vcat(hspool1, hspool2, hspool3, hspool4)
    end
    #
    Zp1 = 0.0
    println(hspool)
    for hscfg in hspool
        bseq = I(2)
        #println(hscfg)
        cfg = [parse(Float64, hscfg[(2*bi-1):2*bi]) for bi in Base.OneTo(Nt)]
        #println(cfg)
        wgt = 1.0
        ηdict = Dict(
            -2 => -3.301360247771569, 
            -1 => -1.049295246550581,
            1 => 1.049295246550581,
            2 => 3.301360247771569
        )
        γdict = Dict(
            -2 => 0.18350341907227408, 
            -1 => 1.8164965809277258,
            1 => 1.8164965809277258,
            2 => 0.18350341907227408
        )
        for bi in Base.OneTo(Nt)
            c = cfg[bi]
            η = ηdict[c]
            γ = γdict[c]
            sqrtma = sqrt(complex(-dt*U/2, 0.0))
            M11 = exp(sqrtma*η)
            M12 = 0.0
            M21 = 0.0
            M22 = exp(sqrtma*η)
            #
            cfg1 = [M11 M12; M21 M22]
            tk = [exp(-dt*mu1) 0; 0 exp(-dt*mu2)]
            bseq = cfg1*tk*bseq
            wgt = wgt*γ*exp(sqrtma*η*(phi-1))
        #println(cfg1)
        #Zp1 += 0.5*det(I(2) + cfg1)
        end
        Zp1 += (0.25^Nt)*det(I(2) + bseq)*wgt
    end
    return Zp1
end


Z1 = Zonesite_tr(4)
println(Z1)

Z2 = Zonesite_hs(4)
println(Z2)

Z3 = Zonesite_hs_Gauge1(4)
@show Z3

Z4 = Zonesite_hs_Quad2(4, 0.5)
println(Z4)
# (phi-1)^2
println(Z4*exp((0.5-1)*(0.5-1)*0.4*1.5/2))
println(bmat_Quad1ADX([0 1; 1 0], Val(-1), 1.5, 0.0, 0.1))
println(bmat_Quad1ADX([0 1; 1 0], Val(2), 1.5, 0.0, 0.1))
println(Δmat_Quad1ND(1, 1, 2, -1, 1.5, 0.0, 0.1))


println(bmat_Quad2ADX([0 1; 1 0], Val(-1), 1.5, 0.0, 0.1))
println(bmat_Quad2ADX([0 1; 1 0], Val(2), 1.5, 0.0, 0.1))
println(Δmat_Quad2ND(1, 1, 2, -1, 1.5, 0.0, 0.1))
mrai = bmat_Quad2ADX([0 1; 1 0], Val(-1), 1.5, 0.0, 0.1)
mat = mrai[1:2, :] + im*mrai[3:4, :]
println((Δmat_Quad2ND(1, 1, 2, -1, 1.5, 0.0, 0.1)[2] + I(2)) * mat)

#fd = fakedata(4)

# + 4exp(-0.1*0.25)
#Z = exp(-0.1*U*0) + 2exp(-0.1*U*0) + exp(-0.1*U*1)


#println(lam)
#Zp1 = 0.5*(1+exp(-lam))*(1+exp(+lam))
#Zp1 = 2*Zp1
##Zp1 = 0.5*(1+exp(-lam-0.05*U))*(1+exp(+lam-0.05*U))
##Zp1 = 2*Zp1
#
#
#
##dup, ddn = fakedata(1)
##Zp1 = 0.5*(1+sum(dup))*(1+sum(ddn))
##Zp1 = 2*Zp1
#
#println(Z)
#println(Zp1)